%% This code runs Stage III of the method.
% assuming a cost shift
close all; clear;
clc;

% define saved files from running this code.
estFile_save = 'est2_output2_c0c1_2007.mat';
fig1_save = 'fig4a_supply_mean2.png';
t1 = 2007; % cost changes from c0 to c1 in year 2007

% estFile_save = 'est2_output2_c0c1_1996.mat';
% fig1_save = 'fig4b_supply_mean2.png';
% t1 = 1996; % cost changes from c0 to c1 in year 1996

%% load data and estimates from Stage I & II
load 'ATE_data_age10.mat'
load 'est1_output_2015.mat' % output year: 1950-2015

Tend = 2015;
ttmax = length(E);
H = (1e+3)*ATE_catch(year<=Tend); % harvest (kg)
Price = R16price(year<=Tend); % price ($/kg)

%% Estimate coefficients in the dynamic Open Access difference eq
fprintf('DeltaE = gamma*(PH) + (gamma*(c0+c1*D))(-E)\n');
% c(t) = c0+c1*D1, cost shift from c0 to c1
D1 = (year >= t1);

% generate variables for regression y = X*b
y = E(2:end)-E(1:end-1); % DeltaE
z = Price(1:end-1).*H(1:end-1); % Price*Catch
x0 = -E(1:end-1);
x1 = -D1(1:end-1).*E(1:end-1);

X = [z x0 x1];
fprintf('y = DeltaE, x1 = PH, x2 = (-E), x1 = (year>=t1)*(-E)\n');
result = fitlm(X,y,'Intercept',false);
disp(result)
[n, p] = size(X);
AIC = n*log(result.SSE/n)+2*p;
fprintf('AIC = %.2f\n',AIC);

b2 = result.Coefficients.Estimate; % point estimates
se2 = result.Coefficients.SE; % std. err
gamma = b2(1);

%% save residuals (for appendix Fig E.1)
% resid2 = y - result.Fitted;
% save('est2_output2_residual_1996.mat','resid2')
%% Monte Carlo 
M1 = 1000; % number of random draws
M2 = 500; % number of MC iterations

rng(54321) % for replication

% random sampling of gamma from Gamma dist.
pd1 = makedist('Gamma','a',(b2(1)/se2(1))^2, 'b',se2(1)^2/b2(1));
beta1_mc = random(pd1,M1,M2); 
% random sampling of gamma*c0 from Gamma dist.
pd2 = makedist('Gamma','a',(b2(2)/se2(2))^2, 'b',se2(2)^2/b2(2));
beta2_mc = random(pd2,M1,M2);
% random sampling of gamma*c1 from Gamma dist.
pd3 = makedist('Gamma','a',(b2(3)/se2(3))^2, 'b',se2(3)^2/b2(3));
beta3_mc = random(pd3,M1,M2);
% sample mean of cost parameter
c0_sample = mean(beta2_mc./beta1_mc);
c1_sample = mean(beta3_mc./beta1_mc);

c0 = mean(c0_sample); se_c0 = std(c0_sample);
c1 = mean(c1_sample); se_c1 = std(c1_sample);

%% save cost and gamma estimates
save(estFile_save,'gamma','c0','se_c0','c1','se_c1')

%% Simulate backward-bending supply using c0 and c1
cost0 = c0; Rscale = R1_scale;
[H_sim0,P_sim0,Hmax0,backP0] = sim_supply(Rscale,q0,Theta,cost0);

cost1 = c1; Rscale = R1_scale;
[H_sim1,P_sim1,Hmax1,backP1] = sim_supply(Rscale,q0,Theta,cost1);

fig1 = figure('Position',[0 0 900 500]);
hold on
h0 = plot((1e-3)*H_sim0,P_sim0,'k','LineWidth',3,...
    'DisplayName','before cost shift');
line([0 (1e-3)*Hmax0],[backP0 backP0],'Color','k','LineStyle',':','LineWidth',2)
h1 = plot((1e-3)*H_sim1,P_sim1,'b-.','LineWidth',3,...
    'DisplayName','after cost shift');
line([0 (1e-3)*Hmax1],[backP1 backP1],'Color','b','LineStyle',':','LineWidth',2)
lgd = legend([h0 h1],'Location','NW','FontSize',12);
lgd.EdgeColor = 'none';
xlabel 'Catch Quantity (metric tons)'
ylabel 'Price ($/kg)'
ylim([0 100])
xlim([0 Inf])
title(sprintf('simulate supply before and after cost shift in %d',t1))
% gtext(sprintf('$%.2f/kg',backP0),'FontSize',12)
% gtext(sprintf('$%.2f/kg',backP1),'FontSize',12,'Color','b')

saveas(fig1,fig1_save)